disp('tjalff');
clear;
clc;

F = @(x, T) -T*log(1+exp(-x/T));
F2 = @(u, m) log((1 + exp(-u)*fn_inveta(m, h0))...
    ./(1 + fn_inveta(m, h0)));

Num_point = 400;
mumu = 1.0;
T = 0.01;
H = 0.001;
% H1 = H - 1e-5;
% H2 = H + 1e-5;
% Nc = 6;
% vpressure = zeros(1, Nc);
% vpressure_dH = zeros(1, Nc);

% n = 3;
% [vk1, vep1, ck1, cX1, cu1] = fn_solve_tba(T,H1,mumu,n,Num_point);
% [vk2, vep2, ck2, cX2, cu2] = fn_solve_tba(T,H2,mumu,n,Num_point);
% [vk, vep, ck, cX, cu] = fn_solve_tba(T,H,mumu,n,Num_point);
% [vk_dH, vep_dH, ck_dH, cu_dH] ...
%     = fn_solve_tba_dH(T, H, vk, vep, ck, cu);
% 
% close(figure(1));
% figure(1);
% plot(vk_dH, vep_dH, '-ro');
% % tvk = [vk1, vk2];
% % tvk = sort(tvk);
% % tvk = unique(tvk);
% tvk = linspace(0, max(vk1(end), vk2(end)), 1000);
% tep1 = interp1(vk1, vep1, tvk, 'pchip');
% tep2 = interp1(vk2, vep2, tvk, 'pchip');
% figure(1); hold on;
% plot(tvk, (tep1 - tep2)/(H1 - H2), '-*');


for n = 6
    
    [vk, vep, ck, cX, cu] = fn_solve_tba(T,H,mumu,n,Num_point);
    
    [vk_dH, vep_dH, ck_dH, cu_dH] ...
    = fn_solve_tba_dH(T, H, vk, vep, ck, cu);
    
    figure(9); hold on;
    plot(fn_ex1(vk), fn_ex2(vep), '-o');
    title('the ep');
    pause(0.01);

    figure(99); hold on;
    plot(fn_ex1(vk_dH), fn_ex2(vep_dH), '-o');
    title('the ep_dH');
    pause(0.01);
    
    fprintf('finished Nc = %d \n', n);
    fprintf('=====================\n\n');
end

cP = cell(1, n);
cQ = cell(1, n);

for ni = 1:n
    
    Pn = fn_u2Pn(ni, H, T, cu{ni});
    cP{ni} = Pn;
    Qn = fn_u2Qn(ni, H, T, cu{ni});
    cQ{ni} = Qn;
    
end

figure(888); hold on;
for ni = 1:n
    
   figure(888); hold on; 
   plot(fn_ex1(ck{ni}), fn_ex2(cX{ni}), '-o'); 
   title('X');
   xlim([-20,20]);
   figure(889); hold on; 
   plot(fn_ex1(ck{ni}), fn_ex2(cP{ni}), '-o');
   title('the P');
   xlim([-20,20]);
   figure(890); hold on;
   plot(fn_ex1(ck{ni}), fn_ex2(cQ{ni}), '-o');
   title('the Q');
   xlim([-20,20]);
    
   figure(220); hold on;
   plot(fn_ex1(ck_dH{ni}), fn_ex2(cu_dH{ni}), '-o');
   xlim([-20,20]);
   
end




% close(figure(1));
% close(figure(2));
% close(figure(3));
% close(figure(4));
% close(figure(5));

% close(figure(1));
% figure(11);
% plot(1:Nc, vpressure, '-o');
% 
% figure(12);
% plot(1:Nc, vpressure_dH, '-o');
% 
% 
% for n = 1:Nc
%     figure(888); hold on;
% %     a = fn_inveta(n, H/T);    
% %     tv = -log((1 + exp(-cu{n})*a)/(1+a));
%     plot(fn_ex1(ck{n}), fn_ex2(cX{n}), '-o');
%     xlim([-40, 40]);
%     title('the X');
%     figure(666); hold on;
%     plot(fn_ex1(ck_dH{n}), fn_ex2(cu_dH{n}), '-o');
%     xlim([-40, 40]);
%     title('the cu_dH');
% end





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%


disp('tjsjlftj');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [vk_dH, vep_dH, ck_dH, cu_dH] ...
    = fn_solve_tba_dH(T, H, vk, vep, ck, cu)

Num_point = 400;

Nc = length(ck);
cP = cell(1, Nc);
cQ = cell(1, Nc);

for n = 1:Nc
    
    Pn = fn_u2Pn(n, H, T, cu{n});
    cP{n} = Pn;
    Qn = fn_u2Qn(n, H, T, cu{n});
    cQ{n} = Qn;
    
end

ck_dH = cell(1, Nc);
cu_dH = cell(1, Nc);

for n = 1:Nc
    
    ck_dH{n} = linspace(0, 1000, 200);
    cu_dH{n} = zeros(1, length(ck_dH{n}));
    
end

vk_dH = vk;
vep_dH = -1/2*tanh(H/2/T)*ones(1, length(vk));

while 1
    
    vep_dH0 = vep_dH;
    ck_dH0 = ck_dH;
    cu_dH0 = cu_dH;
    
    tvep_dH = zeros(1, length(vk));
    for i = 1:length(vk_dH)
        
        k = vk_dH(i);
        res = -1/2*tanh(H/2/T);
%         res = -1;
        res = res + fn_int_an(2, vk, fn_G(vep, T).*vep_dH, k);
            
        for n = 1:Nc
            
            [tx2, tf2] = fn_combine(ck{n}, cP{n},  cQ{n}, ck_dH{n}, cu_dH{n});
            res = res - T*fn_int_an(n, tx2, tf2, k);
            
        end
        
        tvep_dH(i) = res;
        
    end
    
    vep_dH = tvep_dH;
    
    vF = -1/T*fn_G(vep, T).*vep_dH;
    
    vnorm = zeros(1, Nc);
    for n = 1:Nc
        
        [vp, vf] ...
            = fn_iteration2(n, vk, vF, ck, cP, cQ, ck_dH, cu_dH, Num_point);
        ck_dH{n} = vp;
        cu_dH{n} = vf;
        
        vnorm(n) = fn_norm_2(ck_dH{n},cu_dH{n},ck_dH0{n},cu_dH0{n});
        
    end
    
    fprintf('diss = %1.2e   ', fn_norm(vep_dH - vep_dH0));
    fprintf('diss = %1.2e   \n', fn_norm(vnorm));
    
    if fn_norm(vep_dH - vep_dH0)<1e-8 && fn_norm(vnorm)<1e-8
        break;
    end
    
end

end


function [vx, vF] = fn_combine(vk, P, Q, vu, u_dH)

[tvk2, tv2] = fn_prod(vk, Q, vu, u_dH);
vk2 = [tvk2, vk];
vx = unique(sort(vk2));
vF = interp1(vk, P, vx, 'pchip') ...
    - interp1(tvk2, tv2, vx, 'pchip');

end


function [vp, vr] ...
    = fn_iteration2(n, vk, vF, ck, cP, cQ, ck_dH, cu_dH, Num_point)

f = @(x) fn_F2(n, x,vk,vF,ck,cP,cQ,ck_dH,cu_dH);

p_left = 0;
p_right = 1000;
vp = linspace(p_left, p_right, 100);
vr = zeros(1, length(vp));

for i = 1:length(vp)
    
    p = vp(i);
    r = f(p);
    vr(i) = r;
    
end

while 1
    
    vd = abs(vr(2:end) - vr(1:end-1));
    [~, max_region] = max(vd);
    ind_left = max_region;
    ind_right = max_region + 1;
    
    p_left = vp(ind_left);
    p_right = vp(ind_right);
    p_mid = 1/2*(p_left + p_right);
    
    new_r = f(p_mid);
    vp = [vp(1:ind_left), p_mid, vp(ind_right:end)];
    vr = [vr(1:ind_left), new_r, vr(ind_right:end)];
    
    if length(vp)>=Num_point
        break;
    end
    
end

end


function res = fn_F2(n, k, vk, vFT, ck, cP, cQ, ck_dH, cu_dH)

Num = length(cP);
res = fn_int_an(n, vk, vFT, k);

for m = 1:Num
    
    [tk, tf] =  fn_combine(ck{m}, cP{m}, cQ{m}, ck_dH{m}, cu_dH{m});
    res = res + fn_int_Tnm(n, m, tk, tf, k); 
    
end

end


function [vx, vf] = fn_prod(vx1, vf1, vx2, vf2)

vx = [vx1, vx2];
vx = sort(vx);
vx = unique(vx);

vf1_ins = interp1(vx1, vf1, vx, 'pchip');
vf2_ins = interp1(vx2, vf2, vx, 'pchip');
vf = vf1_ins.*vf2_ins;

end


function [vk, vep] = fn_insert(vk0, vep0, T)

Num2 = 200;

x1 = fn_find_zero(vk0, vep0 + 16*T);
% x2 = fn_find_zero(vx, vep + 8*T);
% x0 = fn_find_zero(vk0, vep0);
% x3 = fn_find_zero(vx, vep - 8*T);
x4 = fn_find_zero(vk0, vep0 - 16*T);

vk_insert = linspace(x1, x4, Num2);

vk = [vk0, vk_insert];
vk = sort(vk);

vep = interp1(vk0, vep0, vk);

end


function [vk, vep, ck, cX, cu] = fn_solve_tba(T,H,mumu,Nc,Num_point)

h0 = H/T;
Num = Num_point;

vk = linspace(0, 3*sqrt(mumu + T*log(2*cosh(h0/2))), Num);
% vep_ori = vk.^2 - mumu - T*log(2*cosh(h0/2));
vep_ori = vk.^2 - mumu - T*log(2*cosh(H/2/T));

F = @(x, T) -T*log(1+exp(-x/T));
% F2 = @(u, m) log((1 + exp(-u)*fn_inveta(m, h0))...
%     ./(1 + fn_inveta(m, h0)));
ck = cell(1, Nc);
cX = cell(1, Nc);
cu = cell(1, Nc);

for n = 1:Nc
    ck{n} = linspace(0, 4, Num);
    cX{n} = zeros(1, Num);
    cu{n} = zeros(1, Num);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vep = vep_ori;
while 1
    
    vep0 = vep;
    ck0 = ck;
    cX0 = cX;
    
    vep = vep_ori + fn_an(2, vk, F(vep, T), vk);
    
    for n = 1:Nc
        vep = vep + T*fn_an(n, ck{n}, cX{n}, vk);
    end
    
    vnorm = zeros(1, Nc);
    for n = 1:Nc
        
        vFT = -1/T*F(vep, T);
        [vp, vr, vu] = fn_iteration(n,h0,vk, vFT, ck, cX, Num);
        ck{n} = vp;
        cX{n} = vr;
        cu{n} = vu;
        
        vnorm(n) = fn_norm_2(ck{n},cX{n},ck0{n},cX0{n});
        
    end
    
    fprintf('diss = %1.2e   ', fn_norm(vep - vep0));
    fprintf('diss = %1.2e   \n', fn_norm(vnorm)); 
    
     if fn_norm(vep - vep0)<1e-8 && fn_norm(vnorm)<1e-8
        break;
    end

end

end





function res = fn_norm_2(vx1, vf1, vx2, vf2)

vx = [vx1, vx2];

x1 = vx1(end);
x2 = vx2(end);

xmin = min(x1, x2);
vx(vx>xmin) = [];
vx = sort(vx);

vf1_ins = interp1(vx1, vf1, vx);
vf2_ins = interp1(vx2, vf2, vx);

res = fn_norm(vf1_ins-vf2_ins);

end



% function vX = fn_u2X(m, h0, vu)
% 
% a = fn_inveta(m, h0);
% vX = -log((1 + a*exp(-vu))./(1 + a));
% 
% end


% function vu = fn_X2u(m, h0, vX)
% 
% a = fn_inveta(m, h0);
% vu = -log(1/a*(exp(-vX)-1) + exp(-vX));
% 
% end


function [res, u] = fn_f(n, k, a, vk, vFT, ck, cX)

Num = length(cX);
res = fn_int_an(n, vk, vFT, k);

for m = 1:Num
    
   res = res - fn_int_Tnm(n, m, ck{m}, cX{m}, k); 
    
end


% if res < 0
%     fprintf('lfhwsj \n');
% end

u = res;
res = -log((1 + exp(-res)*a)/(1+a));

end


function [vp, vr, vu] = fn_iteration(n,h0,vk,vFT,ck,cX,Num_point)

a = fn_inveta(n,h0);
f = @(x) fn_f(n, x, a, vk, vFT, ck, cX); 

p_left = 0;
p_right = 1000;

% vp = [p_left, p_right];
% [r_left, u_left] = f(p_left);
% [r_right, u_right] = f(p_right);
% vr = [r_left, r_right];
% vu = [u_left, u_right];

vp = linspace(p_left, p_right, 100);
vr = zeros(1, length(vp));
vu = zeros(1, length(vp));

for i = 1:length(vp)
    p = vp(i);
    [r, u] = f(p);
    vr(i) = r;
    vu(i) = u;
end

while 1
    
    vd = abs(vr(2:end) - vr(1:end-1));
    
    [~, max_region] = max(vd);
    
    ind_left = max_region;
    ind_right = max_region + 1;
    
    p_left = vp(ind_left);
    p_right = vp(ind_right);
    p_mid = 1/2*(p_left + p_right);
    
    [new_r, new_u] = f(p_mid);
    
    vp = [vp(1:ind_left), p_mid, vp(ind_right:end)];
    vr = [vr(1:ind_left), new_r, vr(ind_right:end)];
    vu = [vu(1:ind_left), new_u, vu(ind_right:end)];

    if length(vp) >= Num_point
        break;
    end

end

end


function res = fn_int_Tnm(n, m, vx, vf, k)

res = 0;

if n == m
    
    for i = 2:2:(2*n-2)
        
       res = res + 2*fn_int_an(i, vx, vf, k);  
        
    end
    
    res = res + fn_int_an(2*n, vx, vf, k);
    
else
    
    for i = (abs(m-n)+2):2:(m+n-2)
        
        res = res + 2*fn_int_an(i, vx, vf, k);
        
    end
    
    res = res + fn_int_an(abs(m-n), vx, vf, k);
    res = res + fn_int_an(m+n, vx, vf, k);
    
end

end


function res = fn_int_an(n, vx, vf, k)

res = fn_int_a1(vx/n, vf, k/n);

end


function res = fn_int_a1(vx, vf, k)

vp_tail = vx(end-9:end);
vf_tail = vf(end-9:end);
v2 = 1./vp_tail.^2;
v4 = 1./vp_tail.^4;
kc = vx(end);
[alpha, beta] = fn_fit(v2, v4, vf_tail);

vp0 = linspace(-4, 4, 1000);
vp0 = vp0.^2.*sign(vp0);
vp1 = vp0 + k;
vp2 = vx;
vp3 = vp0 - k;

vp = [vp1, vp2, vp3];
vp(vp<0) = [];
vp(vp>kc) = [];
vp = unique(sort(vp));

vf2 = interp1(vx, vf, vp, 'pchip');
va = 1/2/pi./(1/4 + (k - vp).^2) + 1/2/pi./(1/4 + (k + vp).^2);

vnt = fn_nt(vp);
res = sum(vnt.*va.*vf2);

res = res + alpha*fn_ff2(k, kc) + beta*fn_ff4(k, kc);

end


function res = fn_G(x, T)

if x/T>33
    
    res = 0;
    
else

    res = 1./(1 + exp(x/T));
    
end

end



% function Pn = fn_Pn(H, T, un)
% 
% [Nc, Num] = size(un);
% Pn = zeros(Nc, Num);
% 
% if H>0
% 
%     for n = 1:Nc
%         
%         vu = un(n,:);
%         P = (1 - exp(-vu))/T.*((n+1)*coth(H/2/T*(n+1)) - coth(H/2/T)) ...
%             .*2.*sinh(H/2/T)^2 ...
%             ./(cosh(H/2/T*(2*n+2)) - exp(-vu) + (exp(-vu) - 1).*cosh(H/2/T*2));
%         Pn(n,:) = P;
%         
%     end
% 
% end
% 
% end
% 
% function Qn = fn_Qn(H, T, un)
% 
% [Nc, Num] = size(un);
% Qn = zeros(Nc, Num);
% h0 = H/T;
% 
% for n = 1:Nc
%     vu = un(n,:);
%     eta = fn_eta(n, h0);  
%     Q = 1./(eta*exp(vu) + 1);
%     Qn(n,:) = Q;
% end
% 
% end


function pressure = fn_int_pressure(vx, vep, T)

vx_nt = fn_nt(vx);
pressure = -2*sum(vx_nt.*fn_F(vep,T))/2/pi;

end


function pressure_dH = fn_int_pressure_dH(vx, vep, vep_dH, T)

vx_nt = fn_nt(vx);
pressure_dH = -2*sum(vx_nt.*vep_dH.*fn_G(vep, T))/2/pi;

end


function res = fn_F(x, T)

if -x/T>33
    res = x;
else
    res = -T*log(1 + exp(-x/T));
end

end



% function vu_out = fn_Tnm(n, m, a, vx, vu, vx_out)
% 
% vu_out = 0;
% 
% vf = log((1 + exp(-vu)*a)./(1 + a));
% 
% if m == n
%     
%     for i = 2:2:(2*n-2)
%         
%         vf_out = fn_anf(i, vx, vf, vx_out);
%         vu_out = vu_out + 2*vf_out;
%         
%     end
%     
%     vu_out = vu_out + fn_anf(2*n, vx, vf, vx_out);
%     
% else
%     
%     vu_out = vu_out + fn_anf(abs(m-n), vx, vf, vx_out);
%     vu_out = vu_out + fn_anf(m+n, vx, vf, vx_out);
%     
%     for i = (abs(m-n)+2):2:(m+n-2)
%         
%         vu_out = vu_out + 2*fn_anf(i, vx, vf, vx_out);
%         
%     end
%     
% end
% 
% end


function res = fn_inveta(m, h0)

if h0>2e-3
    res = sinh(h0/2)^2/(sinh(h0/2*(m+1))^2 - sinh(h0/2)^2);
else
    res = 1/(m^2 + 2*m);
end
     
end


% function res = fn_eta(m, h0)
% 
% if h0>0
%     res = (sinh(h0/2*(m+1))^2 - sinh(h0/2)^2)/sinh(h0/2)^2;
% else
%     res = m^2 + 2*m;
% end
% 
% end


function vf_out4 = fn_an(n, vx, vf, vx_out)

vf_out4 = zeros(1, length(vx_out));

for i = 1:length(vx_out)
    
   k = vx_out(i);
   vf_out4(i) = fn_int_an(n, vx, vf, k);
    
end

end



function vf_out4 = fn_anf(n, vx, vf, vx_out)

an = @(x) 1/2/pi*n./((n/2).^2 + x.^2);

vk = vx_out;
vp = vx;

mk = transpose(vk)*ones(1,length(vp));
mp = ones(length(vk),1)*vp;
ma = an(mk - mp) + an(mk + mp);

vpd = vp(2:end) - vp(1:end-1);
vp_nt = zeros(1,length(vp));
vp_nt(1:end-1) = 1/2*vpd;
vp_nt(2:end) = vp_nt(2:end) + 1/2*vpd;

vf_out = ma*diag(vp_nt)*transpose(vf);
vf_out = transpose(vf_out);
vf_out0 = vf_out;

kc = vp(end);
v_tail = vp(end-9:end);
v2 = 1./v_tail.^2;
v4 = 1./v_tail.^4;
vf_tail = vf(end-9:end);
[alpha, beta] = fn_fit(v2, v4, vf_tail);

%  fprintf('alpha = %f, beta = %f \n', alpha, beta);

vf_out_addon2 = alpha ...
    *8/2/pi/n^2*(fn_f2(2*vk/n, 2*kc/n) + fn_f2(-2*vk/n, 2*kc/n));

vf_out_addon4 = beta ...
    *32/2/pi/n^4*(fn_f4(2*vk/n, 2*kc/n) + fn_f4(-2*vk/n, 2*kc/n));

%  vf_out2 = vf_out0 + vf_out_addon2;

vf_out4 = vf_out0 + vf_out_addon2 + vf_out_addon4;

end


function vx_ex1 = fn_ex1(vx)

vx_ex1 = [-vx(end:-1:2), vx];

end

function vx_ex2 = fn_ex2(vx)

vx_ex2 = [vx(end:-1:2), vx];

end


function res = fn_ff2(k, kc)

res = 2/pi/kc/(4*k^2+1) ...
    +4/pi*(4*k^2-1)/(4*k^2+1)^2*(pi/2+atan(2*(k-kc))) ...
    +8*k/pi/(4*k^2+1)^2*log((1+4*(k-kc)^2)/(4*kc^2));

end

function res = fn_ff4(k, kc)

res = 2/3/pi/(1+4*k^2)/kc^3 ...
    +8*k/pi/(1+4*k^2)^2/kc^2 ...
    +8*(48*k^4+8*k^2-1)/pi/(1+4*k^2)^4/kc ...
    +16/pi*(16*k^4-24*k^2+1)/(1+4*k^2)^4*(pi/2+atan(2*(k-kc))) ...
    +64*k/pi*(4*k^2-1)/(4*k^2+1)^4*log((1+4*(k-kc)^2)/(4*kc^2));

end


function res = fn_f2(x, kc)

res = -pi/2*(1 - x.^2)./(1 + x.^2).^2 ...
    + 1/kc./(1 + x.^2) + (1 - x.^2)./(1 + x.^2).^2.*atan(kc - x) ...
    + x./(1 + x.^2).^2.*log((1 + (kc - x).^2)./kc^2);

end


function res = fn_f4(x, kc)

res = 3*pi*(1 - 6*x.^2 + x.^4) ...
    + 6*(1 - 6*x.^2 + x.^4).*atan(x - kc) ...
    + 12*x.*(x.^2 - 1).*log((1 + (x - kc).^2)./kc^2) ...
    + 2*(1 + x.^2).^3/kc^3 ...
    + 6*x.*(1 + x.^2).^2/kc^2 ...
    + 6*(1 + x.^2).*(3*x.^2 - 1)/kc;

res = res./6./(1 + x.^2).^4;

end



function x0 = fn_find_zero(vx, vep)

[a0, a2, a4] = fn_fit_poly(vx, vep);
x0 = sqrt((sqrt(a2^2-4*a0*a4)-a2)/2/a4);

end



function [a0, a2, a4] = fn_fit_poly(vx, vep)

n_right = find(vep>0, 1);
n_left = n_right - 1;

vind = [(n_left - 4):n_left, n_right:(n_right + 4)];

x = vx(vind);
y = vep(vind);

Num = length(x);

V0 = ones(1,Num);
V2 = x.^2;
V4 = x.^4;

n = @(Vec1, Vec2) Vec1*transpose(Vec2);

M = [n(V0,V0), n(V0,V2), n(V0,V4); ...
    n(V2,V0), n(V2,V2), n(V2,V4); ...
    n(V4,V0), n(V4,V2), n(V4,V4)];
V = [n(V0,y); n(V2,y); n(V4,y)];

tv = M\V;
a0 = tv(1);
a2 = tv(2);
a4 = tv(3);

end


function [alpha, beta] = fn_fit(v2, v4, v)

m = [v2*transpose(v2), v2*transpose(v4); ...
    v4*transpose(v2), v4*transpose(v4)];

RHS = [v2*transpose(v); v4*transpose(v)];

tv = m\RHS;
alpha = tv(1);
beta = tv(2);

end


function res = fn_norm(vx)

res = mean(abs(reshape(vx, [1, numel(vx)])));

end


function vx_nt = fn_nt(vx)

vx_nt = zeros(1,length(vx));
vx_len = vx(2:end) - vx(1:end-1);
vx_nt(1:end-1) = 1/2*vx_len;
vx_nt(2:end) = vx_nt(2:end) + 1/2*vx_len;

end


function Pn = fn_u2Pn(n, H, T, un)

if abs(H/T)>2e-3

Pn = (1-exp(-un))/T ...
    .*2*((n+1)*coth(H/2/T*(n+1)) - coth(H/2/T)).*sinh(H/2/T)^2 ...
    ./(cosh(H/T*(n+1)) - exp(-un) - (1-exp(-un))*cosh(H/T));

else
    
    Pn = (1-exp(-un))*(n^2+2*n)/6/T*H/T./(n^2+2*n+exp(-un));
    
end

end


function Qn = fn_u2Qn(n, H, T, un)

a = fn_inveta(n, H/T);
Qn = a./(exp(un) + a);

end

